www.gusucode.com > 超声波测量以及形成图像 对相关信号进行模拟仿真 > 超声波测量以及形成图像 对相关信号进行模拟仿真/digital holograpy/prog/DBDLC.m
function [ U2 ] = DBDLC( U1,z,dxy1,t,center,RPN,varargin ) %DBDLC Diffraction By Discrete Linear Convolution % % Syntax: % U2 = DBDLC(U1,z,dxy1,t,center,RPN,'PropertyName','PropertyValue',... ); % % U1 : the wavefront on the x1 plane % U2 : a part of the free space diffraction of U1 on x2 plane % size of U1 and U2 must be even % z : the distance between these two planes % dxy1 : the sampling interval of U1 % dxy2 : the sampling interval of U2 % t : the parameter to decide dxy2, dxy2=dxy1/t, t must be an integer % center : 1 by 2 array, it is the position of U2's center % center=[centerx,centery] % RPN : 1 by 2 array, it is the point number of U2 % (U2 is a RPN(1) by RPN(2) array) % RPN=[RPNy,RPNx] <== ATTENTION % %-------------------------------------------------------------------------- % PropertyName and PropertyValue: % % inclination {k} | rs | off % it decides whether and which inclination factor is considered in the % integrand % k - use Kirchhoff inclination factor % rs - use Rayleigh-Sommerfeld inclination factor % off - inclination factor is not considered % % interpolation {zero} | sinc % the interpolation method of U1 % zero - interpolate with zeros, in this case, the program is equal % to direct integral of Kirchhoff diffraction formula % sinc - sinc interpolation by inverse fourier transform of an array % which is formed by pading the frequency of U1 with zeros % ------------------------------------------------------------------------ error(nargchk(6,10,nargin)) for n=1:length(varargin) if ~ischar(varargin{n}) error('Property names and values must be characters') end end parray=[0,0]; property=struct('name',{'inclination','interpolate'}); property(1).value={'k','rs','off',''}; property(2).value={'zero','sinc',''}; for l=1:2:length(varargin) namefound=0; valuefound=0; m=0; n=0; while ~namefound && m<length(parray) m=m+1; if strcmp(varargin{l},property(m).name) namefound=1; end end while namefound && ~valuefound && n<length(property(m).value) n=n+1; if strcmp(varargin{l+1},property(m).value{n}) valuefound=1; parray(m)=n; end end if namefound==0 error('wrong property name') elseif valuefound==0 error('wrong property value') end end % ------------------------------------------------------------------------ load lambda; k=2*pi/lambda; [M,N]=size(U1); M=M*t; N=N*t; dxy2=dxy1/t; if parray(2)==0 || parray(2)==1 || parray(2)==3 temp=U1; U1=zeros(M,N); U1(1:t:M-t+1,1:t:N-t+1)=temp; clear temp elseif parray(2)==2 FU1=fourier(U1); FU1=paste(zeros(M,N),FU1); U1=invfourier(FU1); end %-------------------------------------------------------------------------- [P,Q]=meshgrid((1-N/2-RPN(2)/2:N/2+RPN(2)/2-1)*dxy2(1)+center(1),(M/2+RPN(1)/2-1:-1:1-M/2-RPN(1)/2)*dxy2(2)+center(2)); r=sqrt(P.^2+Q.^2+z.^2); clear P Q if parray(1)==0 || parray(1)==1 || parray(1)==4 h=exp(i*k*r)./r.*(z./r+1)/2; elseif parray(1)==2 h=exp(i*k*r)*z./r.^2; else h=exp(i*k*r)./r; end clear r H=fft2(h); clear h %-------------------------------------------------------------------------- U1=paste(zeros(M+RPN(1)-1,N+RPN(2)-1),U1,M/2+1,N/2+1); F=fft2(U1); clear U1 %-------------------------------------------------------------------------- F=F.*H; clear H %-------------------------------------------------------------------------- [P,Q]=meshgrid(0:N+RPN(2)-2,0:M+RPN(1)-2); fp=exp(i*2*pi*((M-1)*Q/(M+RPN(1)-1)+(N-1)*P/(N+RPN(2)-1))); clear P Q %-------------------------------------------------------------------------- F=F.*fp; clear fp %-------------------------------------------------------------------------- U2=ifft2(F); U2=U2(1:RPN(1),1:RPN(2)); U2=-i*U2/lambda*dxy1(1)*dxy1(2);